function R=randtn(m,s,a,b,n)
    A=(a-m)./s; B=(b-m)./s; 
    U=rand(1,n); 
    Delta=.5.*(1+erf(B./sqrt(2)))-.5.*(1+erf(A./sqrt(2))); 
    R=s.*sqrt(2).*erfinv(2.*(Delta.*U+.5.*(1+erf(A./sqrt(2))))-1)+m; 
end